Prediction of ligand binding mode among multiple cross-docking poses by molecular dynamics simulations.

We propose a method to identify the correct binding mode of a ligand with a protein among multiple predicted docking poses. Our method consists of two steps. First, five independent MD simulations with different initial velocities are performed for each docking pose, in order to evaluate its stability. If the root-mean-square deviations (RMSDs) of heavy atoms from the docking pose are larger than a given threshold (2.0 Å) in all five parallel runs, the pose is filtered out and discarded. Then, we perform accurate all-atom binding free energy calculations for the residual poses only. The pose with the lowest binding free energy is identified as the correct pose. As a test case, we applied our method to a previously built cross-docking test set, which included 104 complex systems. We found that the present method could successfully identify the correct ligand binding mode for 72% (75/104) of the complexes for current test set. The possible reasons for the failure of the method in the other cases were investigated in detail, to enable future improvements.

Uncovering abnormal changes in logP after fluorination using molecular dynamics simulations.

The fluorination-induced changes in the logP (1-octanol/water partition coefficient) of ligands were examined by molecular dynamics simulations. The protocol and force field parameters were first evaluated by calculating the logP values for n-alkanes, and their monofluorinated and monochlorinated analogs. Then, the logP values of several test sets (1-butanol, 3-propyl-1H-indole, and analogs fluorinated at the terminal methyl group) were calculated. The calculated results agree well with experiment, and the root mean square error values are 0.61 and 0.68 log units for the GAFF and GAFF2 force fields, respectively. Finally, the logP estimation was extended to a drug molecule, TAK-438, for which fluorination-induced abnormal logP reduction has been observed experimentally. This abnormal change was qualitatively reproduced by the molecular dynamics simulations. We found that the abnormal logP reduction can be mainly attributed to the effect of fluorination-induced dipole change. Our results suggest that molecular simulation is a useful strategy to predict the fluorination-induced change in logP for drug discovery applications.

Multiscale enhanced sampling of glucokinase: Regulation of the enzymatic reaction via a large scale domain motion.

Enhanced sampling yields a comprehensive structural ensemble or a free energy landscape, which is beyond the capability of a conventional molecular dynamics simulation. Our recently developed multiscale enhanced sampling (MSES) method employs a coarse-grained model coupled with the target physical system for the efficient acceleration of the dynamics. MSES has demonstrated applicability to large protein systems in solution, such as intrinsically disordered proteins and protein-protein and protein-ligand interactions. Here, we applied the MSES simulation to an important drug discovery target, glucokinase (GCK), to elucidate the structural basis of the positive cooperativity of the enzymatic reaction at an atomistic resolution. MSES enabled us to compare two sets of the free energy landscapes of GCK, for the glucose-bound and glucose-unbound forms, and thus demonstrated the drastic change of the free energy surface depending on the glucose concentration. In the glucose-bound form, we found two distinct basins separated by a high energy barrier originating from the domain motion and the folding/unfolding of the α13 helix. By contrast, in the glucose-unbound form, a single flat basin extended to the open and super-open states. These features illustrated the two distinct phases achieving the cooperativity, the fast reaction cycle staying in the closed state at a high glucose concentration and the slow cycle primarily in the open/super-open state at a low concentration. The weighted ensemble simulations revealed the kinetics of the structural changes in GCK with the synergetic use of the MSES results; the rate constant of the transition between the closed state and the open/super-open states, kC/O = 1.1 ms-1, is on the same order as the experimental catalytic rate, kcat = 0.22 ms-1. Finally, we discuss the pharmacological activities of GCK activators (small molecular drugs modulating the GCK activity) in terms of the slight changes in the domain motion, depending on their chemical structures as regulators. The present study demonstrated the capability of the enhanced sampling and the associated kinetic calculations for understanding the atomistic structural dynamics of protein systems in physiological environments.

Discovery of 3-Benzyl-1-( trans-4-((5-cyanopyridin-2-yl)amino)cyclohexyl)-1-arylurea Derivatives as Novel and Selective Cyclin-Dependent Kinase 12 (CDK12) Inhibitors.

Cyclin-dependent kinase 12 (CDK12) plays a key role in the coordination of transcription with elongation and mRNA processing. CDK12 mutations found in tumors and CDK12 inhibition sensitize cancer cells to DNA-damaging reagents and DNA-repair inhibitors. This suggests that CDK12 inhibitors are potential therapeutics for cancer that may cause synthetic lethality. Here, we report the discovery of 3-benzyl-1-( trans-4-((5-cyanopyridin-2-yl)amino)cyclohexyl)-1-arylurea derivatives as novel and selective CDK12 inhibitors. Structure-activity relationship studies of a HTS hit, structure-based drug design, and conformation-oriented design using the Cambridge Structural Database afforded the optimized compound 2, which exhibited not only potent CDK12 (and CDK13) inhibitory activity and excellent selectivity but also good physicochemical properties. Furthermore, 2 inhibited the phosphorylation of Ser2 in the C-terminal domain of RNA polymerase II and induced growth inhibition in SK-BR-3 cells. Therefore, 2 represents an excellent chemical probe for functional studies of CDK12 and could be a promising lead compound for drug discovery.

Discovery of a Novel Series of Pyrazolo[1,5-a]pyrimidine-Based Phosphodiesterase 2A Inhibitors Structurally Different from N-((1S)-1-(3-Fluoro-4-(trifluoromethoxy)phenyl)-2-methoxyethyl)-7-methoxy-2-oxo-2,3-dihydropyrido[2,3-b]pyrazine-4(1H)-carboxamide (TAK-915), for the Treatment of Cognitive Disorders.

It has been hypothesized that selective inhibition of phosphodiesterase (PDE) 2A could potentially be a novel approach to treat cognitive impairment in neuropsychiatric and neurodegenerative disorders through augmentation of cyclic nucleotide signaling pathways in brain regions associated with learning and memory. Following our earlier work, this article describes a drug design strategy for a new series of lead compounds structurally distinct from our clinical candidate 2 (TAK-915), and subsequent medicinal chemistry efforts to optimize potency, selectivity over other PDE families, and other preclinical properties including in vitro phototoxicity and in vivo rat plasma clearance. These efforts resulted in the discovery of N-((1S)-2-hydroxy-2-methyl-1-(4-(trifluoromethoxy)phenyl)propyl)-6-methyl-5-(3-methyl-1H-1,2,4-triazol-1-yl)pyrazolo[1,5-a]pyrimidine-3-carboxamide (20), which robustly increased 3',5'-cyclic guanosine monophosphate (cGMP) levels in the rat brain following an oral dose, and moreover, attenuated MK-801-induced episodic memory deficits in a passive avoidance task in rats. These data provide further support to the potential therapeutic utility of PDE2A inhibitors in enhancing cognitive performance.

Exploring the Stability of Ligand Binding Modes to Proteins by Molecular Dynamics Simulations: A Cross-docking Study.

Docking has become an indispensable approach in drug discovery research to predict the binding mode of a ligand. One great challenge in docking is to efficiently refine the correct pose from various putative docking poses through scoring functions. We recently examined the stability of self-docking poses under molecular dynamics (MD) simulations and showed that equilibrium MD simulations have some capability to discriminate between correct and decoy poses. Here, we have extended our previous work to cross-docking studies for practical applications. Three target proteins (thrombin, heat shock protein 90-alpha, and cyclin-dependent kinase 2) of pharmaceutical interest were selected. Three comparable poses (one correct pose and two decoys) for each ligand were then selected from the docking poses. To obtain the docking poses for the three target proteins, we used three different protocols, namely: normal docking, induced fit docking (IFD), and IFD against the homology model. Finally, five parallel MD equilibrium runs were performed on each pose for the statistical analysis. The results showed that the correct poses were generally more stable than the decoy poses under MD. The discrimination capability of MD depends on the strategy. The safest way was to judge a pose as being stable if any one run among five parallel runs was stable under MD. In this case, 95% of the correct poses were retained under MD, and about 25-44% of the decoys could be excluded by the simulations for all cases. On the other hand, if we judge a pose as being stable when any two or three runs were stable, with the risk of incorrectly excluding some correct poses, approximately 31-53% or 39-56% of the two decoys could be excluded by MD, respectively. Our results suggest that simple equilibrium simulations can serve as an effective filter to exclude decoy poses that cannot be distinguished by docking scores from the computationally expensive free-energy calculations.

Discovery of Clinical Candidate N-((1S)-1-(3-Fluoro-4-(trifluoromethoxy)phenyl)-2-methoxyethyl)-7-methoxy-2-oxo-2,3-dihydropyrido[2,3-b]pyrazine-4(1H)-carboxamide (TAK-915): A Highly Potent, Selective, and Brain-Penetrating Phosphodiesterase 2A Inhibitor for the Treatment of Cognitive Disorders.

Phosphodiesterase (PDE) 2A inhibitors have emerged as a novel mechanism with potential therapeutic option to ameliorate cognitive dysfunction in schizophrenia or Alzheimer's disease through upregulation of cyclic nucleotides in the brain and thereby achieve potentiation of cyclic nucleotide signaling pathways. This article details the expedited optimization of our recently disclosed pyrazolo[1,5-a]pyrimidine lead compound 4b, leading to the discovery of clinical candidate 36 (TAK-915), which demonstrates an appropriate combination of potency, PDE selectivity, and favorable pharmacokinetic (PK) properties, including brain penetration. Successful identification of 36 was realized through application of structure-based drug design (SBDD) to further improve potency and PDE selectivity, coupled with prospective design focused on physicochemical properties to deliver brain penetration. Oral administration of 36 demonstrated significant elevation of 3',5'-cyclic guanosine monophosphate (cGMP) levels in mouse brains and improved cognitive performance in a novel object recognition task in rats. Consequently, compound 36 was advanced into human clinical trials.

Exploring the stability of ligand binding modes to proteins by molecular dynamics simulations.

The binding mode prediction is of great importance to structure-based drug design. The discrimination of various binding poses of ligand generated by docking is a great challenge not only to docking score functions but also to the relatively expensive free energy calculation methods. Here we systematically analyzed the stability of various ligand poses under molecular dynamics (MD) simulation. First, a data set of 120 complexes was built based on the typical physicochemical properties of drug-like ligands. Three potential binding poses (one correct pose and two decoys) were selected for each ligand from self-docking in addition to the experimental pose. Then, five independent MD simulations for each pose were performed with different initial velocities for the statistical analysis. Finally, the stabilities of ligand poses under MD were evaluated and compared with the native one from crystal structure. We found that about 94% of the native poses were maintained stable during the simulations, which suggests that MD simulations are accurate enough to judge most experimental binding poses as stable properly. Interestingly, incorrect decoy poses were maintained much less and 38-44% of decoys could be excluded just by performing equilibrium MD simulations, though 56-62% of decoys were stable. The computationally-heavy binding free energy calculation can be performed only for these survived poses.

Design and synthesis of potent and selective pyridazin-4(1H)-one-based PDE10A inhibitors interacting with Tyr683 in the PDE10A selectivity pocket.

Utilizing structure-based drug design techniques, we designed and synthesized phosphodiesterase 10A (PDE10A) inhibitors based on pyridazin-4(1H)-one. These compounds can interact with Tyr683 in the PDE10A selectivity pocket. Pyridazin-4(1H)-one derivative 1 was linked with a benzimidazole group through an alkyl spacer to interact with the OH of Tyr683 and fill the PDE10A selectivity pocket. After optimizing the linker length, we identified 1-(cyclopropylmethyl)-5-[3-(1-methyl-1H-benzimidazol-2-yl)propoxy]-3-(1-phenyl-1H-pyrazol-5-yl)pyridazin-4(1H)-one (16f) as having highly potent PDE10A inhibitory activity (IC50=0.76nM) and perfect selectivity against other PDEs (>13,000-fold, IC50=>10,000nM). The crystal structure of 16f bound to PDE10A revealed that the benzimidazole moiety was located deep within the PDE10A selectivity pocket and interacted with Tyr683. Additionally, a bidentate interaction existed between the 5-alkoxypyridazin-4(1H)-one moiety and the conserved Gln716 present in all PDEs.

Design and synthesis of a novel 2-oxindole scaffold as a highly potent and brain-penetrant phosphodiesterase 10A inhibitor.

Highly potent and brain-penetrant phosphodiesterase 10A (PDE10A) inhibitors based on the 2-oxindole scaffold were designed and synthesized. (2-Oxo-1,3-oxazolidin-3-yl)phenyl derivative 1 showed the high P-glycoprotein (P-gp) efflux (efflux ratio (ER)=6.2) despite the potent PDE10A inhibitory activity (IC50=0.94 nM). We performed an optimization study to improve both the P-gp efflux ratio and PDE10A inhibitory activity by utilizing structure-based drug design (SBDD) techniques based on the X-ray crystal structure with PDE10A. Finally, 1-(cyclopropylmethyl)-4-fluoro-5-[5-methoxy-4-oxo-3-(1-phenyl-1H-pyrazol-5-yl)pyridazin-1(4H)-yl]-3,3-dimethyl-1,3-dihydro-2H-indol-2-one (19e) was identified with improved P-gp efflux (ER=1.4) and an excellent PDE10A inhibitory activity (IC50=0.080 nM). Compound 19e also exhibited satisfactory brain penetration, and suppressed PCP-induced hyperlocomotion with a minimum effective dose of 0.3mg/kg by oral administration in mice.

Discovery of 1-[2-fluoro-4-(1H-pyrazol-1-yl)phenyl]-5-methoxy-3-(1-phenyl-1H-pyrazol-5-yl)pyridazin-4(1H)-one (TAK-063), a highly potent, selective, and orally active phosphodiesterase 10A (PDE10A) inhibitor.

A novel series of pyridazinone-based phosphodiesterase 10A (PDE10A) inhibitors were synthesized. Our optimization efforts using structure-based drug design (SBDD) techniques on the basis of the X-ray crystal structure of PDE10A in complex with hit compound 1 (IC50 = 23 nM; 110-fold selectivity over other PDEs) led to the identification of 1-[2-fluoro-4-(1H-pyrazol-1-yl)phenyl]-5-methoxy-3-(1-phenyl-1H-pyrazol-5-yl)pyridazin-4(1H)-one (27h). Compound 27h has potent inhibitory activity (IC50 = 0.30 nM), excellent selectivity (>15000-fold selectivity over other PDEs), and favorable pharmacokinetics, including high brain penetration, in mice. Oral administration of compound 27h to mice elevated striatal 3',5'-cyclic adenosine monophosphate (cAMP) and 3',5'-cyclic guanosine monophosphate (cGMP) levels at 0.3 mg/kg and showed potent suppression of phencyclidine (PCP)-induced hyperlocomotion at a minimum effective dose (MED) of 0.3 mg/kg. Compound 27h (TAK-063) is currently being evaluated in clinical trials for the treatment of schizophrenia.

Prediction of Ligand Binding Affinity by the Combination of Replica-Exchange Method and Double-Decoupling Method.

A prediction method for ligand binding affinities to proteins is proposed. We first predict the structures of protein-ligand complex by the replica-exchange umbrella sampling or its extension. We then calculate ligand binding affinities based on these predicted ligand-protein bound structures by the double-decoupling method. As a test of the effectiveness of the proposed method, we applied it to the system of the oncoprotein MDM2 and a ligand. The value of the predicted binding affinity turned out to be in good agreement with that from experiments.

Solvation free energies of alanine peptides: the effect of flexibility.

The electrostatic (ΔGel), van der Waals cavity-formation (ΔGvdw), and total (ΔG) solvation free energies for 10 alanine peptides ranging in length (n) from 1 to 10 monomers were calculated. The free energies were computed both with fixed, extended conformations of the peptides and again for some of the peptides without constraints. The solvation free energies, ΔGel, and components ΔGvdw, and ΔG, were found to be linear in n, with the slopes of the best-fit lines being γel, γvdw, and γ, respectively. Both γel and γ were negative for fixed and flexible peptides, and γvdw was negative for fixed peptides. That γvdw was negative was surprising, as experimental data on alkanes, theoretical models, and MD computations on small molecules and model systems generally suggest that γvdw should be positive. A negative γvdw seemingly contradicts the notion that ΔGvdw drives the initial collapse of the protein when it folds by favoring conformations with small surface areas. When we computed ΔGvdw for the flexible peptides, thereby allowing the peptides to assume natural ensembles of more compact conformations, γvdw was positive. Because most proteins do not assume extended conformations, a ΔGvdw that increases with increasing surface area may be typical for globular proteins. An alternative hypothesis is that the collapse is driven by intramolecular interactions. We find few intramolecular H-bonds but show that the intramolecular van der Waals interaction energy is more favorable for the flexible than for the extended peptides, seemingly favoring this hypothesis. The large fluctuations in the vdw energy may make attributing the collapse of the peptide to this intramolecular energy difficult.

Two-dimensional replica-exchange method for predicting protein-ligand binding structures.

We have developed a two-dimensional replica-exchange method for the prediction of protein-ligand binding structures. The first dimension is the umbrella sampling along the reaction coordinate, which is the distance between a protein binding pocket and a ligand. The second dimension is the solute tempering, in which the interaction between a ligand and a protein and water is weakened. The second dimension is introduced to make a ligand follow the umbrella potential more easily and enhance the binding events, which should improve the sampling efficiency. As test cases, we applied our method to two protein-ligand complex systems (MDM2 and HSP 90-alpha). Starting from the configuration in which the protein and the ligand are far away from each other in each system, our method predicted the ligand binding structures in excellent agreement with the experimental data from Protein Data Bank much faster with the improved sampling efficiency than the replica-exchange umbrella sampling method that we have previously developed.

Ligand docking simulations by generalized-ensemble algorithms.

In protein chemistry and structural biology, conventional simulations in physical statistical mechanical ensembles, such as the canonical ensemble with fixed temperature and isobaric-isothermal ensemble with fixed temperature and pressure, face a great difficulty. This is because there exist a huge number of local-minimum-energy states in the system and the conventional simulations tend to get trapped in these states, giving wrong results. Generalized-ensemble algorithms are based on artificial unphysical ensembles and overcome the above difficulty by performing random walks in potential energy, volume, and other physical quantities or their corresponding conjugate parameters such as temperature and pressure. The advantage of generalized-ensemble simulations lies in the fact that they not only avoid getting trapped in states of energy local minima but also allow the calculations of physical quantities as functions of temperature or other parameters from a single simulation run. In this chapter, we review the generalized-ensemble algorithms. Some of their specific examples such as replica-exchange molecular dynamics and replica-exchange umbrella sampling are described in detail. Examples of their applications to drug design are presented.

Prediction of Protein-Ligand Binding Structures by Replica-Exchange Umbrella Sampling Simulations: Application to Kinase Systems.

We have applied our prediction method, which is based on the replica-exchange umbrella sampling for protein-ligand binding structures, to two kinase systems (p38 and JNK3) with two different ligand molecules for each kinase. Starting from configurations in which the protein and the ligand are far away from each other, our method predicted the ligand binding structures in excellent agreement with the experimental data from PDB in all four cases, which suggests the general applicability of our method to kinase systems. In addition, the protein flexibility was shown to be essential to predict the correct binding structure for one of the systems, where dihydroquinolinone was bound to p38 alpha kinase (PDB ID: 1OVE ).

Fast Calculations of Electrostatic Solvation Free Energy from Reconstructed Solvent Density using proximal Radial Distribution Functions.

Although detailed atomic models may be applied for a full description of solvation, simpler phenomenological models are particularly useful to interpret the results for scanning many, large, complex systems where a full atomic model is too computationally expensive to use. Among the most costly are solvation free energy evaluations by simulation. Here we develop a fast way to calculate electrostatic solvation free energy while retaining much of the accuracy of explicit solvent free energy simulation. The basis of our method is to treat the solvent not as a structureless dielectric continuum, but as a structured medium by making use of universal proximal radial distribution functions. Using a deca-alanine peptide as a test case, we compare the use of our theory with free energy simulations and traditional continuum estimates of the electrostatic solvation free energy.

Ab initio prediction of protein-ligand binding structures by replica-exchange umbrella sampling simulations.

We have developed a prediction method for the binding structures of ligands with proteins. Our method consists of three steps. First, replica-exchange umbrella sampling simulations are performed along the distance between a putative binding site of a protein and a ligand as the reaction coordinate. Second, we obtain the potential of mean force (PMF) of the unbiased system using the weighted histogram analysis method and determine the distance that corresponds to the global minimum of PMF. Third, structures that have this global-minimum distance and energy values around the average potential energy are collected and analyzed using the principal component analysis. We predict the binding structure as the global-minimum free energy state on the free energy landscapes along the two major principal component axes. As test cases, we applied our method to five protein-ligand complex systems. Starting from the configuration in which the protein and the ligand are far away from each other in each system, our method predicted the ligand binding structures in excellent agreement with the experimental data from Protein Data Bank.

Peptide conformational preferences in osmolyte solutions: transfer free energies of decaalanine.

The nature in which the protecting osmolyte trimethylamine N-oxide (TMAO) and the denaturing osmolyte urea affect protein stability is investigated, simulating a decaalanine peptide model in multiple conformations of the denatured ensemble. Binary solutions of both osmolytes and mixed osmolyte solutions at physiologically relevant concentrations of 2:1 (urea:TMAO) are studied using standard molecular dynamics simulations and solvation free energy calculations. Component analysis reveals the differences in the importance of the van der Waals (vdW) and electrostatic interactions for protecting and denaturing osmolytes. We find that urea denaturation governed by transfer free energy differences is dominated by vdW attractions, whereas TMAO exerts its effect by causing unfavorable electrostatic interactions both in the binary solution and mixed osmolyte solution. Analysis of the results showed no evidence in the ternary solution of disruption of the correlations among the peptide and osmolytes, nor of significant changes in the strength of the water hydrogen bond network.